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I. INTRODUCTION 

Recent development in laser technology has opened up the possibility of exploring previously 
unattainable regimes of laser-plasma interactions, [|IJ. Intensities of the order of 10 18 W/cm 2 
and higher can now be achieved, implying that goals like compact sources for x-ray lasers, 
0, the fast ignitor concept for inertial confinement fusion, (ICF) || and laser-plasma based 
accelerators, |4j], might soon be within reach. However, a major effort is still required both 
numerically and analytically in order to understand the nonlinear phenomena which arise in 
the presence of such extremely high intensities of the electromagnetic radiation. Besides, a 
good analytical insight is needed in order make numerical simulations possible and to interpret 
their results, ]5]. 

One of the problems which have received particular attention is the combined effect of 
relativistic and striction (ponderomotive) nonlinearities, which occur in the propagation of 
superintense laser pulses through underdense plasmas (i.e., plasmas with uj p < cu, where u p = 
(47m e e 2 /m e ) 1 / 2 is the plasma frequency and lu is the laser carrier frequency), ||. This problem 
is not fully understood yet and there is need of a self-consistent analytical description which 
does not violate the global charge conservation and plasma quasineutrality when describing 
self- focusing and self-channeling, M. 
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A common feature in the above mentioned advanced schemes is that transport of laser ra- 
diation over considerable distances, well beyond the diffraction limit, and without significant 
energy losses is required. In achieving this goal, nonlinear self-focusing and self-channeling play 
an important role. Under the action of an intense laser pulse, electrons tend to be redistributed 
in the transverse direction as an effect of the ponderomotive pressure, the self- channeling phe- 
nomenon. The subsequent self-modification of the radial profile of the refractive index is at the 
origin of the nonlinear self-focusing and filamentation of the laser pulse. Nonlinear self-focusing 
and self-channeling offer a possibility for optical guiding of laser pulses in underdense pulses 
as, for instance, the underdense corona of an ICF target, through the formation of "hollow 
channels" , |J , as experimentally observed by several groups, [j5|,|TU| . 



As was shown in [IT], relativistic self-focusing shows qualitatively new features for ultra- 
intense lasers. The ponderomotive force of super-strong fields expels electrons thus producing 
"vacuum channels" which guide the radiation and stable channeling with power higher than 
the critical one can take place, [I2]]. However, as noticed in [13||, analytical descriptions led to 



the appearance of negative electron densities. This problem was solved by setting the channel 
boundary positions exactly at the point where the electron density became zero. Feit et a!., 
0, showed that this procedure did not conserve the global charge and proposed to include the 
electron temperature effect, which however was not self-consistently evaluated (the tempera- 
ture was assumed to be derived from experimental conditions). Recently, we noticed that for 
an overdense plasma a self-consistent description of self-induced transparency is possible which 
automatically takes into account global charge conservation through Poisson's equation, flHf|. 



The strong analogies between ID overdense and 2D underdense plasmas, allow for an exact 
analysis of the stationary stage of electron cavitation due to the joint effects of relativistic 
and striction nonlinearities in underdense plasmas. This analysis leads to an exact analytical 
description of the transverse structures generated by relativistic self-focusing and also demon- 
strates its multi-filament nature. 

Depending on the incident power and intensity distribution, several qualitatively different solu- 
tions may occur. Our aim is to give an exact analytical description of the stationary stage of the 



fundamental configurations. We will show that, if the incident power is relatively low and the 
intensity has its maximum on-axis, the plasma will react by generating a stationary stage with 
a single channel acting as an optical guide for the propagating radiation. If the incident power 
is increased and the intensity instead has its minimum on-axis (a higher order laser mode), then 
the final stationary stage will display two symmetric channels. Finally, for even higher incident 
powers and maximum on-axis incident intensities, three channels will be generated and so on, 
with the critical power for channel formation depending on the unperturbed plasma density 
and the wave number of the propagating radiation. 

Of particular interest will be the fact that these structures can be interpreted as the final sta- 
tionary stage of the filamentation instability, as shown by the numerical simulations presented 
in §. 

In this paper, we introduce, in Section |TJ, the model equations and the approximations we will 
use to describe electron cavitation in a two-dimensional underdense plasma. The commonly 
presented analysis is discussed in Section [Hi] and our results for single and multi-channel struc- 
tures in 2D, planar geometry are presented in the following Sections, following a brief discussion 
of the physical mechanisms behind the generation of such structures. Finally, some conclusions 



are summarized in Section VI 



II. THE MODEL 

Let us consider the propagation of electromagnetic radiation in a homogeneous plasma. A 
complete description is based on Maxwell's equations for the propagating laser radiation plus a 
model describing the plasma response. For the problem of interest, an important simplification 
comes from the physical context. Considering short pulses with length r such that io~l < r C 



LO, 



pi 



, the ion dynamics can be neglected, ||12|| . Furthermore, we will not be concerned with wake 



field generation, |13| , since the pulse is long enough to allow us to disregard longitudinal charge 
separation, 0. 

Finally, all thermal effects will be disregarded since, at these high intensities, electrons are 
driven to relativistic velocities in just a few optical cycles and the electron pressure gradient is 
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negligible compared to the ponderomotive pressure, [|T5|1 . 

These assumptions define the model we are using to describe our plasma: The ions are consid- 
ered as an immobile neutralizing background and the electrons as a cold relativistic fluid. Our 
set of self-consistent equations derived from Maxwell's equations and the equation of motion 
for the electron component, assuming the Coulomb gauge, reads 

^ 2a 1 d 2 A 1<9 4tt 

V 2 A-^ = - S V„ + -AW (1) 

VV = 4vre(iV - iV ) (2) 

r/ryv = -A + Vip (3) 
c 

^ = e<p - mc 2 (>y - 1) (4) 

V • A = 0. (5) 

Here 7 = l/^/l — v 2 /c 2 is the relativistic factor, N is the electron density, Nq is the equilibrium 
density, — e and m the electron charge and mass respectively, A is the electromagnetic vector 
potential, (p is the electrostatic scalar potential and ip is a scalar function which expresses the 
electron canonical momentum. Details on the derivation of this model can be found in 
Eqs.(||) and (f|) imply that we are assuming vortex- free motion of the electrons. Taking the 
divergence of Eq.([l|) and using Eqs.(g) and (f|) we find that the charge conservation law 

dN 

— + V-(iVv) = (6) 

is automatically satisfied, i.e., the total charge is conserved. However, when dealing with 
necessarily simplified models describing the stationary regime in presence of electron cavitation 
phenomena, the condition of plasma quasineutrality is not obviously conserved, fTl"Hi"5|. This 
point must be carefully discussed when constructing new solutions and it will lead to the 
breaking of the Hamiltonian model, thus allowing for multi-filament structures. 

The assumption made on the pulse duration implies that the electron fluid has time to 
approach a quasi-steady state, (Tj|. Therefore, it is interesting to describe what kind of sta- 
tionary state the system will reach, neglecting any transient phenomena. This leads to further 
simplification, because we can neglect the time dependence in Eq.(|4]) and, since v before the 
passage of the laser pulse must be zero, it follows from Eqs.(|3|) and (f|) that tp = 0. 
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In order to single out the fast optical time scale we adopt the slowly varying envelope 
approximation, factorizing the normalized vector potential as: 

eA 



mc 2 



a±(rj_)r_i_ exp [i(hz — cut)] + c.c. (7) 



Assuming the paraxial approximation k± << h, where k± is the transverse component of the 
laser wavenumber and h is the propagation constant, the parallel component of the vector 
potential is negligible if compared to the transverse ones and the incident radiation may be 
assumed to be circularly polarized without loss of generality In what follows, we will drop 
the subscript denoting the perpendicular component of the various quantities. The resulting 
system of equations, after a few algebraic manipulations, is: 

V 2 a+ fl-^)a = (8) 

V 2 = a(n - 1) (9) 

= 7 — 1 if and only if n ^ (10) 

7 = Vl + a 2 (11) 

where 

° = T^m (12) 

k = Lo/c is the vacuum wavenumber and we have introduced the following normalization: 
n = N /N cr with N cr = mu 2 / (4ne 2 ), n = N/N , <p = etp/(mc 2 ), r = k^Jl - h 2 /k 2 f±. 

III. COMMON ANALYSIS 

Let us consider a two dimensional geometry for a plasma extending in the z direction, i.e. 
along the laser propagation direction. We will restrict ourselves to a ID transverse model in 



order to emphasize the main features of multi-filament structures. In Eqs.(§)-(|TT|) the prop- 
agation constant h plays the role of a free parameter which, together with the background 
plasma density n a , defines what kind of filament structures can be realized as a final state of 
the self-focusing evolution. In reality, it would depend on several parameters and factors as 



the laser power, the geometrical configuration (the angle of focusing, for instance) and on the 
prehistory of the process. 



The self-channeling we are interested in is realized only when a > 1, i.e., for underdense 
plasmas when 1 > h/k > y/1 — no and for overdense plasmas with no > 1 when h/k < 1. 

The complete mathematical analogy between the present model and the one introduced in 
17] suggests that our plasma will react to the laser action with the formation of regions depleted 



of electrons, where the laser electromagnetic radiation is trapped, a consequence of the well 
known phenomenon of electron cavitation and channeling. Electrons tend to be expelled from 
the focal spot by the laser ponderomotive force and, at the same time under such extreme 
conditions, they acquire relativistic quiver velocities. These effects both contribute to a self- 
induced modification of the radial profile of the refractive index and a consequent nonlinear 
trapping of the laser radiation in finite plasma regions. This modification is the basic mechanism 
in the optical guiding of laser pulses in plasmas. It is possible to give an exact analytical 
description of the asymptotic stationary plasma-field structures generated in the transverse 
direction for different values of a. As we will see, these structures consist of one or more 
channels, depending on the corresponding incident power. The most delicate point in the 
analysis, will be the determination of such structures complying with global charge conservation. 



Eqs. (|8|)-(|TT|) were analyzed in detail in [[L3|] and also in [[TJ[] with respect to both underdense and 
overdense plasmas and solutions were found in the form of continuous functions. Fundamental 
to those analysis is the Hamiltonian structure of the set of equations (|8|)-(|TT|), which reads 
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He = 2(1 + a 2 / ' ~ 2 {2aVTT ^ ~ ° 2) (13) 
where the prime denotes the derivative with respect to the transverse coordinate x. As n(x) — > 1 
and both a(x) and a'(x) vanish for x — * oo, the integral of motion equals 

He = H E0 = -a. (14) 

It follows that there is an exact soliton-like analytical solution given by 

A m cosh [|e | 1/2 (^ - x (0) )l 

a ( x ) = : —TTnZ 17^ r-7 ( 15 ) 
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where €q = 1 — a and the parameter defines the peak position of the function ([15]) which is 
given by A m = 2 [a (a — l)] 1 ^ 2 . Once this solution is known, we also have a description for the 
electron density through Poisson's equation (|9|) and the equation of motion (|IUD: 

n = 3(1 + a 2 ) + 2 V1 + Q2 (H e - a 2 ) (16) 

The minimum electron density in a cavity is given by n m i n = 1 — 4(cv — l) 2 , which implies that, 
for a > 1.5, this solution leads to the unphysical result of a negative electron density. Therefore, 
if a < 1.5, i.e., for propagation constants lying in the interval ^1 — 2no/3 < h/k > y/1 — n , 
we have solutions of expressed by the continuous functions (|i5|) and (|16|). 
The field and density structures and the corresponding power, related to the propagation con- 
stant, are presented in Fig. |l] and |2|. 



a=1.4 
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FIG. 1. Plasma-field structures (dashed line) and electron density distribution (continuous line) for 
the case of a single filament in an under dense plasma, for a fixed value of a less then 1.5. In this case 
a = 1.4. 




a 

FIG. 2. The power needed to generate a single filament structure is shown as a function of the 
parameter a for a less then 1.5. All quantities are dimensionless. 
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It should be emphasized that, for a < 1.5, since the system is fully described by the 
Hamiltonian (13), there are no other structures except this single-filament ones. 

The important question is what will happen for higher incident powers or, in other words, 



for a > 1.5? The procedure for constructing a solution followed in [|T|,|T3|,|T^] , which consisted 
in assuming the electron density to vanish within the interval where the solution for the density 



(16) is negative, led to non conservation of the global charge. However, what is happening is 
that the ponderomotive force is pushing electrons away from the central axis, while the force 
due to charge separation acts in the opposite direction. Thus, when an equilibrium is reached, 
we have the formation of a stationary structure consisting of a channel emptied of its electrons. 
This means that the global structure of the solution consists of two parts, the first one described 
by the Hamiltonian (|H|) while the second, describing the depletion regions where the electron 
density vanishes, has the typical vacuum Hamiltonian: 

1 



(17) 



In Fig. |3] the phase portrait of the full system is presented for the single filament case with 
a = 2 while in Fig. H the same phase portrait is shown for a more complicated multi-filament 



case, with a 



a=2. 
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FIG. 3. Phase portrait for a system that develops a single-filament structure. The thick line 
represents the trajectory followed by the system, starting from the separatrix at a = 0, a' = 0, then 
moving on along the vacuum trajectory to finally come back to the starting point along the separatrix 
again which represents the symmetrical plasma region. The corresponding plasma-field structures are 
illustrated in Figj5|. 

a= 2. 




FIG. 4. Phase portrait for a more complicated case with multiple filaments. The thick line 
represents the trajectory in the case of four filaments. Solid lines represent trajectories relative to 
electron layers while dashed lines are relative to depletion regions and the large dashed line separates 



regions with positive and negative electron density, as follows from Eq. (16). 



The continuous soliton-like solution described by ( |I5D corresponds here to the separatrix 
trajectory and its starting and final point is a = 0, a' = 0. As pointed out, this solution breaks 
down for higher values of a and we indicate on the phase portrait the curve beyond which the 
electron density (|i~6l) formally becomes negative. 

Beyond the limit curve for the electron density we have to introduce the "vacuum" part of the 
solution. Our system has left the separatrix and has started to move along the vacuum one. 
The boundary position up to which the electrons are displaced is determined by the equilibrium 
condition between the two forces acting on them, as described by the equation of motion 

<P' = i (is) 



and by the conservation of the total charge, which means that, in order to conserve the total 
charge, the boundary positions can now be determined by inserting the equilibrium condition 
([[ID into Poisson's equation and integrating it. 



IV. SINGLE-FILAMENT SOLUTIONS 

Let us consider a localized solution with one peak for the intensity. Its structure is defined 
by the closed trajectory (0 — 1 — 1' — 0) shown in Fig. 0(a). We will treat this kind of solution 
as a single-filament solution. 

Integrating Poisson's equation over the whole interval we obtain 

x d = ---^= (19) 

where is the field amplitude at the boundary position Xd- At the same time, we have to 
match the field in the vacuum channel 

a(x) = Ay cosx (20) 

and its first derivative to the field and its first derivative in the plasma region, that is, 

A 2 v = a d + a' d 2 (21) 

and 



V 



x d = - arctan I — (22) 

Given the integral of motion He = from the two equations for Xd we obtain a transcendental 
equation for the boundary amplitude ad 



tan ( v - }=— n [ 2c *W 1 + a d ~ !) - a d\ 1 (23) 

OC J (Xd 

which can be solved numerically, so that now we know everything about the structures generated 
in this case, see Fig. ||. It is important to be careful when solving this equation since it 
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has multiple solutions but we have to choose the first which satisfies the condition of charge 
conservation. 

ar= 2.4 
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FIG. 5. Plasma- field structures (dashed line) and electron density distribution (continuous line) for 
the case of a single channel for a fixed value of a = 2. All quantities are dimensionless. 



The calculation of the total power for this single channel configuration as a function of a 
is straightforward: P = a 2 (x)dx, and the result is presented in FigJ^. For comparison, 
we present here also the total power calculated within the previous model, when the boundary 
position Xd was assumed to be the one where the electron density vanished. In this case, 
the total charge being not conserved, there was an excess of positive charge which led to a 
much higher power required in order to overcome the restoring force due to this charge excess. 
Consequently, the power needed to generate such structures was overestimated. Besides, it 
is interesting to see how, for increasing values of a and consequently increasing values of the 
required power, the width of the central vacuum channel becomes larger but, after the initial 
rapid growth, its increase becomes slower, see the box in Fig.0. 
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FIG. 6. Total power (continuous line) and channel width w = 2xd (box) for the single filament case. 
For comparison, the dashed line shows the total power calculated according to the commonly used 
model, without taking into account global charge conservation. 



V. MULTI-FILAMENT STRUCTURES 



It is evident from our analysis that, due to the requirements of global charge conservation 
and to the symmetry imposed with respect to the z-axis, for a fixed a the single filament 
configuration and the power necessary to generate this structure are uniquely determined. If 
this power is exceeded, the incident electromagnetic radiation is strong enough to spread along 
the transverse direction of the plasma channel over a distance larger than in the previous 
case, but still finite. For a sufficiently strong power, the final stationary state will present a 
multiple channel structure which can again be analytically described. The power necessary for 
generating each of these structures can be exactly calculated as well. As we will show, this 
allows for the definition of a threshold power for the generation of multi-filament structures. 

Let us start with the case of an intensity distribution with a minimum on the symmetry 
axis. An example of the trajectory in the phase space for the double channel case is given in 
Fig. 0(b), indicated as (0 — 2 — 2' — 3' — 3 — 0). To construct the field structure we can start 
from inside the plasma region at x — > +oo, where we know the integral of motion and the 
expression for the decaying field and the electron density, see Eqs. flT5| ) and flH3D respectively. 
When we come to the depletion region there is a certain freedom in the choice of the boundary 
amplitude a^, as we are going backwards from the last plasma region towards the central axis. 
The only requirement for ad is that the electron density must not be negative, therefore we can 
fix the boundary amplitude (and therefore the boundary position x rf as well) to any value up 
to a maximum for which the density at the boundary vanishes. This means that, for a fixed 
a the two peak solution is not unique and there is a certain power range for generating such a 
structure. For a given the field in the vacuum region, 



is completely determined from the matching conditions, but now the vacuum channel extends 




(24) 
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from the Xd to a certain x\ (see Fig. ^) which is to be determined taking into account charge 
conservation. 

In order to construct a structure with only one degree of complexity more than for the single 
channel, we stop at the next plasma layer, which will be centered on the symmetry axis. An 
analytical expression for the field in this central plasma layer can be derived by solving the 
equation for the vector potential. Now the solution is not localized as before, therefore the 
boundary conditions and the Hamiltonian He = Hex > are not known. The solution is 
expressed in terms of two-parameter elliptic functions as 

2qcn[e\ /2 (x-xW)] 



2+[( 9 2+i)i/2_i] sn 2[ e V2( a ._ a .(i) )] > El - ~, ^ 

2qsn[((e 1 +l) 2 -a 2 ) 1 / 2 (x-xW)/2] 
q 2 — sn 2 



a < Tim < a, 



a(x\ 



i%ei+i) a -a 2 )V2 (x-xW)/2] ' > a 



where e x = (a 2 + 1 + 2H E i) 1/2 , q = [(ei + a) 2 - 1] 1/2 and q = + a + l)/{e 1 + a- l)] 1 / 2 , 
while k = [{a 2 - ( £l - l) 2 )/4e 1 } 1 / 2 and k = [{ £l - l) 2 - a 2 )/{e x + l) 2 - a 2 )] 1 ' 2 are the moduli 
of the elliptic integrals of the first kind respectively for the two cases. These solutions were 
presented in [H] for the problem of self-induced transparency of an overdense plasma. Imposing 
the conservation of the total charge by integrating Poisson's equation from x = to x = oo 
with the equilibrium condition defined by the equation of motion, see Eq. flTB|) , we obtain a 
transcendental equation for the quantity £ = X\ — x^- 

£ = 9(£)-g(o) (26) 

where 

m = y«+6.)i = _ arotan K) (27) 

2a [1 + Af / cos 2 (£ + £ )] 1/2 \adj 
The solution of this equation gives a complete description of the vacuum layer since x,i is 
already determined. A necessary condition for this equation to have a non-trivial solution is 
that cf{£ = -£o) > 1, i.e. 

9 '<« = -*>> = saw 5 >h (28) 

which cannot be satisfied unless a > 1.5. This leads to the conclusion that if a < 1.5, that is, if 
the propagation constant of the wave vector h is not large enough, our system will never reach 
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a stationary state and it will display only a dynamical behaviour with the electromagnetic per- 
turbation propagating along the transverse direction. Otherwise, we can numerically calculate 
how the minimum boundary intensity, such that Eq. ( p6|) still has a solution, depends on a and 
also find a solution of this equation and calculate x\. To finally obtain a complete description 
of the central plasma layer we only have to apply the boundary conditions at x\ to determine 
the parameters that are still unknown, TCei and while the symmetry axis is determined 
as the symmetry axis of the elliptic function. In Fig.|^ a double-channel structure is shown for 
fixed a and for maximum amplitude at the last boundary (so that the electron density at this 
boundary vanishes). It should be noticed that, for a fixed value of the boundary amplitude 
ad, the width of the vacuum channels and the peak intensity in these channels increase with 
a. Furthermore, the maximum possible boundary amplitude itself is an increasing function 
of a and, for any given value of this parameter, such a maximum amplitude determines the 
maximum power we can deliver to the plasma in order to generate a double-channel structure. 
Exceeding this maximum power will force the system to generate a structure with one more 
filament and therefore we can talk about a threshold power for the generation of multi-filament 
structures. 

ar= 2. 
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FIG. 7. Double-channel structure for a plasma with a fixed a = 2 and maximum possible intensity 
a d at the last boundary (continuous line) and electron density distribution (dotted line). All quantities 
are dimensionless. 

In Fig.§ we show Protai calculated for varying a and for fixed to its maximum possible 
value. In the box it is shown instead how the total power varies with the intensity at the last 
boundary. 
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FIG. 8. Total power required to generate a double channel structure versus a for the case of max- 
imum possible intensity at the last boundary and (box) as a function of the boundary intensity for 
fixed a = 2. 



The apparently anomalous behaviour for low boundary intensities is due to the fact that 
the left branch corresponds to a new kind of two-filament solution whose phase portrait field 
structures are presented in Fig. 0. 



As can be seen following the trajectory (0 — 1 — 2 — 3 — 2/ — 1/ — 0), in this case the 
field amplitude, once it leaves the separatrix, never crosses the zero point until it reaches the 
separatrix again. For those periodic trajectories lying inside the separatrix, the Hamiltonian 
TCei is less than —a and the field solution for the electron layers assume a different form and 
are now described in terms of elliptic functions as 

( , = (g - fi ~ ~ (el /2 /e 2 )s^[e 2 (x - X W)]V /2 (9Q) 
{) l-2(e 1 /a + e 1 -l)sn 2 [£ 2 (x-x( 1 ))] { } 

where E\ = (a 2 + 1 + Hsoo)°°^ e , £2 = [« 2 — (ei — l) 2 ]/2 and the module of the elliptic integral 

1 /2 

of the first kind is k = e 1 /e 2 . The procedure to define the electron cavitation boundaries is 
the same as the one followed previously to build the structures presented in Fig.[7|. 
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FIG. 9. (a) Phase portrait for a system that develops a double-filament structure. In this case 
the field amplitude never vanishes apart from ±oo on the separatrix. The corresponding plasma-field 
structures are illustrated in part (b). 

The solution we have constructed and the corresponding choice of a closed one-cycle trajec- 
tory in the phase space is not unique. We can pass a depletion region not only at the point 3/ 
to form a one cycle trajectory, but also at the point b' in order to create a periodic trajectory 
as b' — b — a — a', see Figs. |3](b) and §(a). Following this trajectory means that we will have 
a new structure with new channels and plasma layers, which were not present in the double 
channel structure previously described. 



We would like to underline the fact that the periodic trajectory shown in Fig. ^ corresponds 
to a particular configuration as the points b f , a', a, b are related by a complete symmetry. It is 
possible to see what this means by looking at the field structures described by such a trajectory 
(see Fig. 110): The central channels are completely symmetric, at each boundary we have the 



same intensity. It is again the request of global charge conservation which leads, by integrating 
Poisson's equation, to a transcendental equation for the field at the boundary of the new plasma 
layer: 
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tan 



a d [2HEi + 2a Jl + a 



7. 2 ,1 V2' 



ft 



■[2^ 1 + 2ft V /l + a2-a2] 1 / 2 



(30) 



This equation is similar to Eq.(|23|) but now the Hamiltonian value 7^^ = TCei is the one defined 
for the new electron layer. In the case of a single-peak field distribution, this equation had a 
unique solution, consequently we can add to the two-filament distribution whole periodic cycles, 
i.e., even numbers of filaments by jumping to consecutive vacuum trajectories at the point bi 



or repeating the same vacuum trajectories, as seen in Fig.|10| where a six-filament structure is 
presented, which corresponds to the trajectory (0 — 2 — 2' — b' — a! — a — b — 3' — 3 — 0), with 
two cycles along the periodic trajectory b' — a! — a — b. Again, it is important to notice that 
this structure, completely symmetric, is peculiar of a planar geometry. 

a= 2. 
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FIG. 10. Six-channel structure for a plasma with fixed a = 2 and maximum possible intensity at 
the last boundary (continuous line) and electron density distribution (dotted line). All quantities are 
dimensionless. 



We can also add to the single filament configuration an odd number of filaments, considering 
trajectories corresponding to a certain number of cycles plus half a cycle. Consider, for example, 
the trajectory (0 — 2 — 2' — b' — b — 4 — 4' — 0), where the points 4, 4' are symmetrical to 2', 2 
respectively, which corresponds to a three-filament structure. The result is shown, for a fixed 
value of ft in Fig.pd|. 
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a= 2. 




FIG. 11. Three-channel structure for a plasma with fixed a = 2 and maximum possible intensity at 
the last boundary (continuous line) and electron density distribution (dotted line). All quantities are 
dimensionless. 



Finally, in Fig.[L2], we present in the same graphic the calculated maximum powers as func- 
tions of a for three of the different cases we have analyzed, single, double and triple channel. 
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FIG. 12. Total power as a function of a for the three different cases: single channel (continuous 
line), double channel (dashed line), triple channel (dotted line). In each case, the amplitude at the 
last boundary was chosen as the maximum possible one. 

The same procedure that we have described, may also be applied to the case presented in 
Fig.|S] where an integer number of filaments can be added since a full cycle occurs within a 
half-space of the phase portrait. 

Therefore, by using this procedure, we can construct multi-filament solutions that exist only 
for a > 1.5 and represent plasma channels with electron cavitation. They differ from each other 
because of the laser power transported along these channels and, as there is a minimum laser 
power required for exciting such structures, we can define the power thresholds for creating 
non-single filament structures. 
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VI. CONCLUSIONS 



In conclusion, we have presented an exact analysis of self-channeling structures generated 
as a consequence of the relativistic self-focusing due to the interaction of ultra-intense laser 
radiation with an underdense plasma. In this analysis, the plasma quasineutrality condition is 
accurately taken into consideration and this quantitatively affects some results on channeling 
laser power. The essential point of the analysis is that it allows us to prove the multi-filament 
nature of the relativistic self-focusing and to calculate the threshold power for exciting multi- 
filament structures. Such a result is not achievable in media with a local nonlinearity, like the 
Kerr one, because the governing equation has an overall Hamiltonian structure. In the case we 
have analyzed, each electron cavitation channel corresponds to a certain part of the trajectory 
followed by the system in the phase space, each part with its own Hamiltonian value, as shown 
in FigfJ. Concerning the definition of a threshold power, it is interesting to notice how filament 
structures with regions depleted of their electrons can be generated by the interaction if the 
parameter a is greater than 1.5. For a < 1.5 we have only single filament field structures, with 
no depletion regions. As soon as a exceeds 1.5 there are plasma regions that are emptied of 
their electrons and the number of filaments thus generated increases with increasing incident 
power. We can therefore define the maximum power incident on a plasma with a = 1.5 as 
the real threshold power for generating non-single filament structures (see Fig.^(b)). The same 
construction procedure followed for the single-filament solution can then be easily extended to a 
more realistic axisymmetrical case to obtain the analogous of Figs. (5) and (7) , while a real 2D 
transverse approach is needed for the multi-filament structure problem, especially for the case 
presented in Fig. (8) where a number of equal filaments have been added to the fundamental 
structure. 
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